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We present a method for extracting gravitational radiation from a three-dimensional numerical 
relativity simulation and, using the extracted data, to provide outer boundary conditions. The 
method treats dynamical gravitational variables as nonspherical perturbations of Schwarzschild ge- 
ometry. We discuss a code which implements this method and present results of tests which have 
been performed with a three dimensional numerical relativity code. 
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Numerical relativity represents the only currently vi- 
able method for obtaining solutions to Einstein equations 
for highly dynamical and strong field sources of gravita- 
tional radiation. Using these techniques to study coa- 
lescing black hole binaries is the purpose of the multi- 
institutional Binary Black Hole "Grand Challenge" Al- 
liance effort which is presently underway in the United 
States. This effort is also motivated by the prospect of 
observations with the next generation of gravitational 
wave detectors. 

In addition to tremendous demands on computational 
resources, implementing the standard 3 + 1 HSj] formula- 
tion of Einstein theory as a Cauchy problem l4] is compli- 
cated considerably by the necessity of imposing bound- 
ary conditions which maintain numerical accuracy and 
the physical correctness of the solution. Both inner and 
outer boundary conditions have received considerable at- 
tention. Recent efforts on interior boundaries have fo- 
cused on the excision of the interior of the black hole from 
the computational domain (see, for example, |^). This 
paper will concentrate on the problem of outer boundary 
conditions applied at a finite radius around a source of 



gravitational waves. 

Proper boundary conditions on spacelike slices of 
asymptotically flat spacetimes are essential for the ac- 
curate computation of the gravitational waveforms pro- 
duced in the strong field region that represent the obser- 
vationally relevant aspect of the computation. Since it is 
not feasible to simulate on spacelike slices out to arbitrar- 
ily large distances from the source, it is necessary to ex- 
tract gravitational waves comparatively near the strong 
field region and to have boundary conditions that allow 
radiation to pass cleanly off the mesh. If poor outgo- 
ing boundary conditions are imposed, spurious radiation 
is produced which can contaminate the computed grav- 
itational waveform. Additionally, the outer boundary is 
usually close enough to the isolated source that backscat- 
ter of radiation from curvature is significant. This source 
of incoming radiation needs to be built into the outer 
boundary conditions. An approach to the extraction of 
gravitational wave information and the computation of 
outer boundary conditions that exploits the matching of 
the interior numerical solution with an exterior perturba- 
tive solution on spacelike slices has been developed dur- 
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ing the past decade and applied to a number of different 
physical scenarios |^^. Extension of these techniques 
to three dimensional (3D) simulations has been one of 
the efforts of the Alliance. A parallel development is also 
underway in the Alliance that matches interior Cauchy 
solutions to exterior solutions on characteristic hypersur- 
faces 1^. 

In this Letter we report the first successful application 
of the perturbative matching approach to provide outer 
boundary conditions for a 3D numerical relativity code. 
In this context, the perturbative matching method should 
be viewed as a general "module" that can be coupled to 
any nonlinear "interior" 3D code. While the latter solves 
the Einstein equations in the highly dynamical strong 
field region, the former extracts gravitational wave data 
and imposes outer boundary conditions. In the follow- 
ing we give an overview of the method and present test 
results obtained from the evolution of linear waves. 

Nonspherical perturbations of Schwarzschild geome- 
try have been studied for many years in the context of 
black hole perturbations [ p"oUTl]] and as a way to extract 
gauge-invariant information about gravitational radia- 
tion 1^. More recently Schwarzschild perturbation the- 
ory has been found to be useful in studying the late-time 
behavior of the coalescence of compact binaries in a nu- 
merical simulation after the horizon has formed ||l^,|l3|,|| . 

We base our treatment of Schwarzschild perturbation 
theory on a recent hyperbolic formulation of Einstein 
field equations jl^. A principal advantage of this ap- 
proach is that we can easily derive perturbative wave 
equations in terms of the standard 3-1-1 variables, mak- 
ing it straightforward to match exterior perturbative so- 
lutions to interior numerical ones. This method leads to 
spatially gauge-invariant radial wave equations for each 
angular mode |l5| which we have used: to extract grav- 
itational radiation from a 3D numerical relativity simula- 
tion, bj to evolve this information to large radius yielding 
an approximate asymptotic waveform and cj to provide 
outer boundary conditions for such a simulation. One of 
the major advantages of the perturbative method is that 
we have replaced the computationally expensive 3D evo- 
lution of the gravitational waves via the nonlinear Ein- 
stein equations with a set of ID linear equations we can 
integrate to high accuracy and minimal computational 
costs. 

We split the gravitational quantities of interest into 
background and perturbed parts: the three-metric gij = 



hij, the extrinsic curvature Kij 



the 



lapse function N — N + a and the shift vector = 
/?' -I- f % where the tilde denotes background quantities. 
We assume a Schwarzschild background, 

(jijdx'dx^ = N-^dr^ + r^{d9^ + sin^ Odc/)^) , (1) 



(It follows that Kij = = /3'). The perturbed parts have 
arbitrary angular dependence. 

We use this approximation to linearize the hyperbolic 
equations and, in particular, we find that the wave equa- 
tion for Kij reduces to a linear wave equation for Kij in- 
volving also the background lapse We separate the 
angular dependence in this equation by expanding Kij in 
terms of tensor spherical harmonics (ei)^ , . . . , {fn)ij. In 
the notation of lO], we have 



ax{t,r){ei)ij +rb^{t, r)(e2)y + 
N-^a+{t,r){h\j + rh+{t,r){h\, + 



(3) 



where ax, 6x are the odd-parity multipoles, while a+, 
6+, c+, c?+ are the even-parity ones. Note that we have 
suppressed the angular indices {i, m) of (ei)ij, . . . , (A)^ 
and of Ox , . . . , d+ over which there is an implicit sum. 

In odd-parity, we take the wave equation for K^g 
and use the linearized momentum constraints, Vfe(K'^j — 
S\k-'j) = 0, to eliminate the odd-parity amplitude (fox)- 
In even-parity, we use the wave equation for K^r to- 
gether with the wave equation obtained for the trace 
K = K'' ■ — h{t,r)Y^^ and the linearized momentum con- 
straints. In this way we eliminate 6+ and c+ and obtain 
two coupled equations for a+ and h. For each {£, m) 
mode, we therefore have one odd-parity equation 
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and two coupled even-parity equations. 
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(The usual Regge- Wheeler and Zerilli equations can be 
obtained by a more complete analysis p^.) 

We now turn to implementation and tests of this tech- 
nique. Consider a numerical relativity simulation which 



2 




FIG. 1. Location of the difFerent outer boundaries and 
of the extraction 2-sphere for two successive timeslices. The 
dark shaded region shows the spatial domain over which the 
3D nonlinear equation are solved. 

evolves Einstein equations in either the standard 3 + 1 
(ADM) or hyperbolic form on a "interior" 3D grid. (No 
restrictions need be put on the choice of the 3D coordi- 
nate system.) In the tests described here, we have used 
the ADM 3D interior code of the Alliance to evolve 
linear, £ = 2, m = 0, Teukolsky waves ||l^]. The interior 
3D grid uses (topologically) Cartesian coordinates and 
up to (129)'^ gridpoints. In the results shown, the ini- 
tial data consists of a metric formed from the sum of the 
background Cartesian metric and a Gaussian envelope of 
amplitude A — 10~^, width b — 1 and of a vanishing ex- 
trinsic curvature. The ADM evolution is undertaken on 
a cubical grid of extent {x, y, z} — ± 4, with unit lapse 
and zero shift. 

During each time step, the procedure for extracting 
radiation and imposing outer boundary conditions pro- 
ceeds in three steps: (1 ) extraction of the independent 
amplitudes and their time derivatives on a 2-sphere of 
radius inscribed in the grid, (2) evolution of the ra- 
dial wave equations (^-(^ out to large radii using the 
extracted amplitudes to construct inner boundary val- 
ues on the "exterior" ID grid and (3) reconstruction and 
"injection" of Kij and dtKij at specified gridpoints at or 
near the boundary of the interior grid. The last step uses 
the momentum constraints and the inverse of the trans- 
formation employed in step (1 ). In Figure ^ we show 
the schematic location of the different outer boundaries 
and of the extraction 2-sphere on two successive spacelike 
time slices at i = to and t = ti. 

The first step involves interpolating Kij and dtKij 
onto the 2-sphere and evaluating their projections onto 
the appropriate conjugate tensor spherical harmonics to 
eliminate the angular dependence and obtain the cor- 
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FIG. 2. Convergence to the analytic solution of the ex- 
tracted (r = 1) and evolved (r — 8) multipole (a+)2o. The 
amplitude is scaled by to compensate for the radial fall-off. 

responding multipoles. The precise location of the 2- 
sphere within the interior 3D grid depends on the specific 
case under investigation. In general, for the perturbative 
matching to be justified, we require that the extraction 
surface be located in a region where the gravitational 
field is close to Schwarzschild. One must be careful that 
the numerical errors introduced by the interior evolution 
do not dominate the exterior solution. (Further details 
will be presented elsewhere fl^.) 

The second step entails evolving equations (0)-(^). In 
the test presented here, we consider a flat background 
with M = 0. The initial data for the ID exterior grid 
is set consistently with the interior 3D initial data. Dur- 
ing each time integration of the ID exterior grid, the 
extracted multipoles are used as inner boundary values 
and standard Sommcrfeld outgoing wave conditions are 
imposed at the outer boundary (e.g., at r = 30). The 
top diagram of Fig. ^ shows the only relevant multipole 
for this test |]l9| , extracted from the interior 3D grid at a 
2-sphere of radius = 1 pO| ; the lower diagram shows 
its evolved waveform at a radius r = 8. Different curves 
refer to different resolutions of the interior 3D grid and 
show the convergence to the analytic solution. 

The third and final step consists of computing outer 
boundary values for the interior 3D grid, and this step 
can proceed in one of two ways. In the first method, the 
values of Kij reconstructed from the exterior ID data 
are injected as Dirichlct data at all the boundary points 
of the interior 3D grid. (No restrictions are put on the 
shape of the 2-surface where the data are injected.) For 
a code using a hyperbolic formulation, dtKij can also 
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FIG. 3. Convergence to zero of the L2 norm of the error 
of integrated over the 3D outer boundary surface. 



be provided at the boundary. Moreover, the evolution 
equation of the three-metric can be integrated using Kij 
at the boundary so Dirichlet data for gij need not be 
provided there. The second method relates, at the 3D 
outer boundary, the null derivatives of the extrinsic cur- 
vature obtained from the interior grid (Kij) and from 
the perturbative module (kij , since background extrinsic 
curvature is assumed to be zero): 



d_ 

dt 



Kij) -f - (Kij - Kij) = (7) 



This method resembles a Sommerfeld condition but is 
more general since it can be used in regions where the 
radiation is not dominated by the asymptotic outgoing 
behavior. Moreover it takes into account arbitrary angu- 
lar dependence, as well as the effects of a Schwarzschild 
black hole background. Experimentation with the two 
methods has shown that this "perturbative Sommerfeld" 
approach is very accurate and generally more stable than 
the direct injection of Dirichlet data. Figure ^ shows 
convergence to zero of the L2 norm of the difference be- 
tween one component of Kij (namely K^z) and its ana- 
lytic value, integrated over the entire outer boundary of 
the interior 3D grid. 

An important property of any outer boundary imple- 
mentation is that it should allow the outgoing radiation 
to escape freely to infinity. In practice, however, a fi- 
nite discretization always produces a certain amount of 
reflection at the outer boundary and this could drive in- 
stabilities which grow exponentially. We have compared 
the long-term stability obtained with the perturbative 
outer boundary module with that of the best "alterna- 
tive" outer boundary conditions we have tried: namely a 



FIG. 4. Long-term evolution of the L2 norm of the error 
in Kzz at the 3D outer boundary surface. The inset shows 
Kzz at the outer boundary along the a;— axis. The interior 
grid has (33)^ gridpoints and the code is stopped when the 
norm is unity. 



standard outgoing Sommerfeld condition. Fig. ^ shows 
the results of this comparison both on a long time scale 
and on a shorter time scale. (The calculations were per- 
formed using a very coarse resolution.) As shown in the 
inset (first echoes at t 7), the amount of reflection 
produced with the perturbative method is smaller than 
that produced by the Sommerfeld condition. Moreover, 
the use of perturbative boundary conditions delays the 
onset of exponential error growth and allows for a much 
longer evolution (see main figure) . Increasing the interior 
resolution and the number of (^, m) modes used we can 
further prolong the running time. Using (49)^ interior 
gridpoints, we are able to evolve the code for more than 
50 crossing times. In strong contrast, when a standard 
Sommerfeld condition is used, increased interior resolu- 
tion results in a shorter running time [Q. 

This work was supported by the NSF Binary Black 
Hole Grand Challenge Grant Nos. NSF PHY 93-18152, 
NSF PHY 93-10083, ASC 93-18152 (ARPA supple- 
mented). Computations were performed at NPAC (Syra- 
cuse University) and at NCSA (University of Illinois at 
Urbana-Champaign). 



Information about the Binary Black Hole Grand Chal- 
lenge Alliance, its goals and present status of research caii 
be found at the web page: tittp : //www.npac . syr . edu/j 



4 



projects/bh/. 

[2] R. Arnowitt, S. Deser and C. W. Misner, in Gravitation, 

edited by L. Witten, Wiley, New York (1962). 
[3] J. W. York in Sources of Gravitational Radiation, edited 

by L. Smarr, Cambridge Univ. Press, Cambridge (1979). 
[4] Y. Choquet-Bruhat and J. W. York in General Relativity 

and Gravitation I, edited by A. Held, Plenum, New York 

(1980). 

[5] E. Seidel and W.-M. Suen, Phys. Rev. Lett., 69, 1845 
(1992). 

[6] A. M. Abrahams and C. R. Evans, Phys. Rev. D, 37, 317 

(1988); Phys. Rev. D, 42, 2585 (1990). 
[7] A. M. Abrahams et al, Phys. Rev. D 45, 3544 (1992); 

P. Anninos et al, Phys. Rev. Lett. 71, 2851 (1993). 
[8] A. M. Abrahams, S. L. Shapiro and S. A. Teukolsky, 

Phys. Rev. D 51, 4295 (1995); A. M. Abrahams and R. H. 

Price, Phys. Rev. D 53, 1963, Ibidem, 1972 (1996). 
[9] N. T. Bishop, R. Gomez, P. R. Holvorcem, R. A. 

Matzner, P. Papadopoulos and J. Winicour Phys. Rev. 

Lett. 76, 4303 (1996); N. T. Bishop, R. Gomez, L. 

Lehner, J. Winicour Phys. Rev. D. 54, 6153 (1996). 
[10] T. Regge and J. A. Wheeler, Phys. Rev. 108, 1063 (1957); 

F. Zerilh, Phys. Rev. Lett. 24 737 (1970). 
[11] V. Moncrief, Annals of Physics 88, 323 (1974). 
[12] R. H. Price and J. PuUin, Phys. Rev. Lett. 72, 3297 

(1994). 

[13] A. M. Abrahams and G. Cook, Phys. Rev. D , (1994). 

[14] Y. Choquet-Bruhat and J. W. York, C. R. Acad. Sci. 
Paris, 321, 1089 (1995); A. Abrahams, A. Anderson, 
Y. Choquet-Bruhat and J. W. York, Phys. Rev. Lett. 
75, 3377 (1995). 

[15] A. M. Abrahams, A. Anderson, C. R. Evans, C. Lea, 
J. Lenaghan, M. Rupright and J. W. York, In prepara- 
tion. 

[16] The Binary Black Hole Grand Challenge Alliance, In 
preparation. 

[17] S. Teukolsky, Phys. Rev. D 26, 745 (1982). 
[18] A. M. Abrahams, L. RezzoUa and M. Rupright, In prepa- 
ration. 

[19] Other {£, m) multipoles are also computed and used in 
the decomposition; but, given the initial data chosen, 
their amplitudes are several orders of magnitude smaller 

[20] On a flat background spacetime, the position of the ex- 
traction 2-sphere is arbitrary. Similar results have been 
obtained also for = 0.5, 1.5, 2, 2.5, 3, 3.5 



5 



